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ABSTRACT 

Using a 2.5D time-dependent numerical code we recently developed, we solve 
the full compressible Navier-Stokes equations to determine the structure of the 
boundary layer between the white dwarf and the accretion disk in non-magnetic 
cataclysmic variable systems. In this preliminary work, our numerical approach 
does not include radiation. In the energy equation, we either take the dissipation 
function ($) into account or we assume that the energy dissipated by viscous pro- 
cesses is instantly radiated away ($ = 0). For a slowly rotating non- magnetized 
accreting white dwarf, the accretion disk extends all the way to the stellar sur- 
face. There, the matter impacts and spreads towards the poles as new matter 
continuously piles up behind it. We carry out numerical simulations for different 
values of the alpha viscosity parameter (a), corresponding to different mass ac- 
cretion rates. In the high viscosity cases {a = 0.1), the spreading boundary layer 
sets off a gravity wave in the surface matter. The accretion flow moves super- 
sonically over the cusp making it susceptible to the rapid development of gravity 
wave and/or Kelvin- Helmholtz shearing instabilities. This BL is optically thick 
and extends more than 30 degrees to either side of the disk plane after only 3/4 
of a Keplerian rotation period {tK=ids). In the low viscosity cases (a = 0.001) 
, the spreading boundary layer does not set off gravity waves and it is optically 
thin. 

Subject headings: accretion, accretion disks - binaries: close novae, cata- 
clysmic variables — white dwarfs — methods: numerical 
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1. Introduction: Accreting White Dwarfs in Cataclysmic Variables 



Cataclysmic variables (CVs) form an interesting class of short-period close binary sys- 
tems compri sing a hot white dwarf (WD) and a rela tively lower mass red dwarf star filling its 
Roche lobe (ICrawford &: Kraft Ill956l : iKraft Ill962l ). In such systems hydrogen- rich matter 
from the red dwarf exits through the inner Lagrange point (LI) and fiows towards the white 
dwarf. In the absence of strong magnetic fields, the matter forms an accretion disk around 
the white dw arf due to the excess angular momentum originating from the orbita l motion of 



the binary (jPrendergast fc Burbidgelll968l : iFlannery Ill974j : iLubow fc Shu Ill975l ). On-going 
accretion at a low rate (quiescence) is interrupted every few weeks to r nonths by in tense ac- 
cretion (outburst) of days to weeks - a dwarf nova (DN) accretion event (IBath Ill972l ). thereby 
increasing the luminosity of the systems by several magnitudes. A thermal instability in the 
accretion disk is believed to trigger the increase in the mass transfe r rate (m) through th e 
disk and thus an increase in the rate of gravitational energy release (jCannizzo et al.lll988l ). 



CV systems are divi ded in sub-classes according to the duration, occurrence an d am- 
plitude of their outburst Jflack fc la Douslll993l : IWarneil E995I : iRitter fc Ko"lb1ll998h : e.g. 
dwarf nova systems (DNs) are non-magnetic and accrete through a disk, they spend most of 
their time in the quiescent state; nova-like systems (NLs) are disk-systems found mostly in 
the high outburst state; polars are devoid of an accretion disk due to their strong magnetic 
fields (the matter is tunneled through the magnetic field lines onto the poles of the WD); 
and intermediate polars (IPs) have truncated inner disks due to their moderately strong 
magnetic fields. Another CV subclass is the classical nova (or just "nova";), charact erized 



by an episode of unstable thermonuclear burning (the thermonuclear runaway - TNR; [Rose 



(119681 )). All CV system are believed to undergo a classical nova explosion every few thousand 
years or more, when enough hydrogen-rich material accumulated in the envelope to reach 
the ignition point at the electron- degenerate base o f t he envelope - where it is compressed 
under the large gravity of the WD ( Starrfield 1971al Jbl: Istarrfield et al. 1972). In the present 
work we concentrate on the study of the non (or weakly) magnetized DN systems, where the 
accretion disk extends all the way to the surface of the WD. 

Since the white dwarf is the most common end-product of stellar evolution 90% of all 
the stars in the Galaxy have evolved or will evolve into white dwarfs) and the accretion disk 
is the most common universal structure resulting from mass transfer with angular momen- 
tum, and both can be directly observed in cataclysmic variable systems (in the ultraviolet), 
an understanding of the accretion process in cataclysmic variables is the first step in a global 
understanding of accretion in other systems throughout the universe. These include young 
stellar objects, accretion onto neutron stars and black holes, and the most difficult to study, 
active galactic nuclei. Therefore, accreting white dwarfs in cataclysmic variables are the best 
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astronomical laboratories to study the physics of accretion disks. 



1.1. The Accretion Disk & Boundary Layer in One-Dimension 



In the disk, masneto-hydrody namic (MHD) turb ulence (IShakura fc Sunyaevlll973l ) due 
to a magneto-rotational instability fiBalbus et al.lll994l ) dissipates potential energy and trans- 
fers angular mom entum outwar d. As a result, the disk matter slowly spirals inwards towards 
the white dwarf (jPringldll98ll ). The total potential energy of accretion can be released at 
the maximum rate 



GM,m 



mr^Vtj^{r^), 



where G is the gravitational constant, is the mass of the WD, r^ its radius, m is the 
mass accretion rate, and VLk{t^) is the Keplerian angular velocity at one stellar radius r^,. 
This is the maximum amount of energy that can be extracted from the accretion process 
per unit of time and the actual luminosity can be smaller than this (see b elow). In the 
standard disk theory (IShakura &: Sunyaevlll973l : iLynden-Bell &: Pringlelll974j ). the accretion 
disk is axisymmetric, geometrically thin in the vertical dimension (it has a thickness H 
such that H/r « 1), and the energy dissipated by the (turbulent) viscosity is instantly 
radiated locally in the -^z directions. Only half of the accretion luminosity is emitted by 
the disk {Ldisk = Laccf^), since the matter is st ill moving at a nearly Keple rian velocity, 
vk ~ \/GM^,/r^, before it is ultimately accreted ( ILynden-Bell Sz Pringldll974j ). 



The remaining accretion energy must, therefore, be dissipate d in order for the matter 
to be accreted onto the surface of the more slowly rotating WD ( lPringlelll98ll ). The disk 
matter is decelerated in that region where the inner disk reaches the stellar surface: the 
boundary layer (BL). The height of the BL can be comparable to the scale height H of the 
accretion disk. The remaining rotational kinetic energy of the accreting matter dissipated 
in the boundary layer per unit of time (Lbl) is nearly equal to half of the total luminosity 
of the accreting matter (Lace), namely: 



Lbl = (1 - f3^) h 



mr 



disk 



*-{Ql{r,)-Q 



') 



(2) 



where f3 = f2^,/r2;^(r^,), and Q^, is the stellar angular velocity. This is so because the material 
accrete d onto the WD su rface corotates with it and keeps a fraction of the rotational kinetic 
energy. iKluzniak I (119871 ). however, pointed out that part of the BL energy goes into spinning 
up the accreting star through the shear at the stellar equator, and hence, the fraction that 
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can be radiated away is smaller than the one given in eq.(2). Namely, one has 

.2 



2 TflT 



(3) 



Therefore, for a non-rotating WD: Lbl = L^isk = Lacc/'2- However, for a rotating WD 
Lbl < Ldisk, and the total energy radiated from the accretion process {Lbl + Ldisk) is 
actually smaller than Lace as defined by eq.(l). The energy "kept" by the corotating accreted 
material per unit of time is : 



(3'L 



mr 



disk 



(4) 



and the power invested to spin up the star through the shear at the equator (eq.2-eq.3) is: 



[(1 - - (1 - py]Ld^sk = 2/5(1 - /?) W = - P)Lacc. 



(5) 



The above relations are correct to order of rj = 1 — Vm/r^,, where Vm is the radius at which 
the gradient of the angular velo city cha i iges s i gn dQ/dr = (for further detail s of t he 
one-dimensional analysis see also iRegev I (119831 ): iKleyl (Il99l[ ): iPopham fc NarayanI (Il995l )). 
Usually 7] is of the order of H/r or smaller. For a WD star rotating at 10% of the Keplerian 
velocity (i.e. several lOOkm/s) the BL is expected to radiate Lbl = O.SlLdisk, while a 
fraction of O.lSLdisk is invested in spinning up the star. For a WD star rotating at 30% of 
the Keplerian velocity (i.e. ~ lOOOkm/s) the BL is expected to radiate Lbl = OAOLdisk- For 
a star rotating near break-up the angular velocity in the disk does not have an extremum 
and decreases steadily outward, therefore the disk is expected to spin-down the fast rotating 
WD (jPopham fc NarayanI 1 19911 ). In this case there is no boundary layer. 



Since the BL is much smaller than the disk itself and each radiate almost equal amounts 
of energy, the disk emits mainly in the optical and near UV, whereas the much h otter BL 
emits in the far UV and in X-rays. During out burst, when m 
Cannizzo et al.lll988f) . observations (e.g. 



1987 



10~^ - 10-W^/yr (IWarner 



Cordova et a"l~l Jl98oh : iMauch et al. I fll995[ ): 



Mauche I (120041 )) reveal that the BL is optically thick and emits mainly in the FUV and 



ory ( 


^rinele & Savoniie 


( 


197d); 


Reeev 


( 


1983^; 


(1995 


); 


Collins et al 


( 


1998[): 


Obach & Glatzel 


10-12 


-10-^^M^/yt{ 


Warner 


1987 


'), observations 



Jl999hl. During quiescence, when m 



( I2OO3I . I2OO5I )) reveal that the boundary layer is optically thin and emits in the hard X-ray 
band with a temperature of the order of lO^K, a s expected by th e theory (e.g. iPringle fc Savoniie 
(ll979h:lTvlerda1 Jl98lh : [King fc Shavivl Jl984[ ): lshaviv I Jl987h : [Naravan fc Popham I Jl993h : 
Popham Ul999h l However m addition, in the outer region of the BL, where the BL meets 
the disk, the optical thickness becomes lar ger again (~ 1) and that region emits in the FUV 
with a temperature Tgjj reaching ~ lO^K ( iPopham 1119991 ). 
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1.2. Observational Background 



Observationally, it was sliown almost three decades ago t hat non-magneti c CVs do emi t 



some fraction of their luminosity in the X-ray 


bands 


Becker & Marshall 


( 


1981 


); 


Cordova & Mason 


(1983) 


ROSAT All-Sky Survey ( 


Beuermann & Thomas 


1993^ 



e^ 



Cordova et al. I (11981 Il98lh : 



standard disk and boundary layer theories, during quiescence ha rd X-rays (lOkeV and higher 



were observed frorn a small region close to the WD (e.g. Ivan Teeseling et al. I (119961 ): 



Mukai et al. I (Il997l )). whil e during outbu rst this emission is replaced by soft X-ray and 
EUV (see e.g. the review of iMauche I (119961 ) and the references therein). The EUV region of 
the spectrum is very difficult to observe because the absorption cross section of (ISM) neutral 
hydrogen is very high. Because of this, only a few system s have been successfully observed 
in the EUV, such as e.g. VW Hyi jMauche et al. I Jl99lh : for which N{H) ~ 6 x lO^^cm'^, 
Pohdan et al"~l Jl99oh l Assuming that the X-ray and EUV emissions are from the BL, and 
the optical and UV emissions a re from the dis k , many previous studies foun d a very low 
ratio LBi/Ldisk (~ 0.001 - 0.04, iMauche et al~l Jl99lh : iHoare fc Drewl Jl99lhl durin g out- 
burst, and a ratio of Lbl/ Ldisk ~ 0.25 in quiescence (e.g. for VW Hyi, iBelloni et al. I (Il99ll ) 
assuming that the WD contributes 50 percent of the UV flux). The X-ray observations of 
underluminous boundary layers culminated with the i^OS'^T observations of ten cataclysmic 
variables (BA Cam, YZ Cnc, CP Com, VW H y i, WX Hyi, TY PsA, V3885 Sgr, CY UMa, 
CY Vel, IX Vel) by Ivan Teeseling fc Verbunt I (Il994l). in which 8 systems were c aught in 
quiescence and 2 systems were caught in outburst. Ivan Teeseling &: Verbunt I (119941 ) derived 
X-ray fluxes from their observations and UV fluxes from existing WE observations and found 
that ratio of the X-ray Luminosity to the UV+Optical Luminosity is much smaller than one. 



Ma ny proce s ses we r e discu s sed to e xplain the " missing boundary l ayer" ( e.g. iFerland et al. 



Shaviv 


( 


1987 


); 


King 


( 


1997 


); 



(I1995I )). The r nain idea was that the kinetic energy of the BL could also be converted 



intowinds (e.g. iKing &: Shaviv 



(Il978h ). heating (e.g. 



(Il984f)'). W D or b el t rotation (e.g. 



Shaviv fc Starrfield 



(fl987h : 



Regey fc Shara 



vected into the outer stellar envelop (e.g.. Godon (1996a, 1997bl ): [Popham (jl997 )). Some 



K^ippenhahn fc Thomas 
( 1989^: or ma ybe ad 



systems have been observed to h ave a WD rot a ting a t a rather large rotational velocity, of 



the order of ~ 1, OOOkm/sec (e.g. ICheng et al. I (119971 ): iPandel et al. I (120051 ) ). which implies 
(from eq.3) LsL/Ldisk = 1/2. However, the observed X-ray luminosities have been much 
smaller than this and would imply a near-Keplerian rotation rate for many systems, which 
is not the case. 



Recent XMM-Newton observations (jPandel et al. II2005I ). taking into account the con- 
tribution of the disk, WD and BL in the optical, UV and X-ray bands, found no evidence of 
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an underluminous BL in 8 quiescent dwarf novae (OY Car, WW Get, AB Dra, U Gem, VW 
Hyi, T Leo, TY PsA, SU UMa) and the data are consistent with Lbl ~ Ldisk- The main 
difference with previous studies was that the WD contribution was taken into account with 
reahstic temperature s taken from the hterature and could dominate over the disk in the UV. 
Pandel et aL I (120051 ) basically assumed Ldisk = Luv + Lopt — Lwd and Lbl = Lx, and found 
LblI Ldisk ~ 1 for 6 objects among 8. For VW Hyi and U Gem (with Lsi/Ldisk ~ 1/3, 1/2 
respectively) they suggest that the inner disk is truncated at r ~ SRwd to explain the 
discrepancy. 



Go don &: SionI (120051 ) further computed the contribution of all the emitting components 



(WD, disk, BL) in the optical, UV and X-ray based oii the e xistin g multiwaveleng th obser- 
vations of VW Hyi, t h e st andar d disk mo del (jPringldll98ll ) and iPopham I (jl999l )'s model 
of the boundary layer. iPopham I (119991 ) has shown that (in quiescence) the BL emits part 
of it energy in the UV band from that region where the outer edge of the BL meets (and 
radiates energy into) the optically thick inner edge of the disk. iGodon fc SionI (120051 ) sug- 
gested that the second component often observed in the quiescent UV spectr a of DNe (the 
so-call ed accretion belt also detected in the FUSE spectra of VW Hyi — (1 Go don et al 



2004bl )) is the UV emission from the outer BL. Taking the UV contribution of this sec- 
ond component {L2nd) into account iGodon fc SionI (120051 ) showed that the luminosity of the 
BL of VW Hyi in quiescence {Lbl = Lx + L2nd) is as ex pected from th e theo ry (namely 
Lbl = Ldisk = Lopt + Luv — Lwd), therefore agreeing with IPandel et al. I (120051 ) that there 
is no missing boundary layer even for VW Hyi and without disk truncation. 

On the other hand, spectroscopic UV observations of accreting WDs have also reached 



a rather advanced stage. Observations have been carried out for a range of CVs. (iLong et al. 



1993 



2004a 



Froning et al.ll2001l : iSzkodv et al.ll2002l : IWelsh et al.ll2003l : iSion et al.ll2005l : IGodon et al. 



to cite just a few) in quiescence and outburst. The observations suggest that the spec- 
trum is made up of several parts, specifically, the underlying WD, the accretion disk, the 
accretion belt that might form on the surface of the star, the hot spot (where the matter 
inflowing from the LI point hits the outer rim of the accretion disk). The temperature of 
each of these components proves to be a most important diagnostic, though as better data 
come in, we can hope that other hydrodynamical variables also become importai it diagnos- 



tics. The case of VW Hydri is especially important because STIS measurements (ISion et al. 



20051 ) have captured that object during its transition to outburst. The VW Hyi transition to 
outburst caught by STIS was apparently an outside-in outburst. The UV lagged the optical. 
The STIS showed that longer FUV wavelengths changed much faster and manifested the 
flux sooner than shorter wavelengths. It has been found that the temperature of the accret- 
ing star can be raised by as much as 50% during outburst. The elevated WD temperature 
returns to its quiescence value soon after the outburst. The elevated surface temperature 
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of the star may, therefore, be a good diagnostic of the energy dissipated as the accretion 
stream impacts, spreads out and slows down and reaches co-rotation with the star. The rate 
at which the WD cools following the heating during the outburst is not only diagnostic of 
the mass/depth of the heated surface layers but also potentially a good diagnostic of the 
optical depth of the accreted material. A detailed understanding of BL physics may also 
help one understand the so-called UV delay where the rise in the UV emission lags the rise in 
the optical by several hours, indicating that t he UV emission might be more symptomatic of 
the physics of the BL (ILivio fc Pringldll992l ). It is important to note here that our present 
numerical simulations are relevant to the changes in the WD temperature observed on a 
time scale of the order of outburst /quiescence cycle (e.g. VW Hyi), and do n ot apply to 



compressional heating taking place on a ni uch longer evolutionary time scale (ISion I Il995 



Godon fc Sionll200l l2003l : IPiro et al. Il2005[ ) 



1.3. The Two-Dimensional Boundary Layer 



The BL has typically been solved in a r nodel where the averaging of the vertical structure 
which is employed in accretion disk studies (IShakura fc Sunyaevll 19731 ) has been extrapolated 
to the surface. One also makes the additional simplifying assumption that the vertical 
velocity in the boundary layer is zero, just as it is in the outer parts of a thin disk mod el. This 
reduces the description of the BL to a one dimensional problem in the radial direction (jPringle 
198ll : iMeyer fc Meyer-Hofmeisteiill982l : iPopham fc Narayarull995l : ICoUins et al.lll998l ) , where 
one actually solves the "slim disk" equations for the BL region (namely, the one- dimensional 
disk equation s -|- radiation in the radial direction; as first proposed in the pioneering work of 
Regev J 1983h'). Howeve r, this assumption is not necessarily correct sinc e the BL may spread 
out (^Ferland et al.l 119821 ) . The best one- dimensional models of the BL ( Popham fc Narayan 
19951 ) predict a rise in the BL temperature during outburst that is much larger than the 
observed one. 



Recognizing the multi-dimensionality of th e problem several au t hors a ttacked the prob 



lem directly using numerical methods (e.g. [Robertson fc FrankI (119861 )) culminating in 



(Kiev & Henslei 


1987; 


Kley 


1989. 


1991) 


Kley 


( 


1989. 


199] 


) were very bold and 



solve the problem has until now precluded following the problem numerically on a long evo- 
lution time scale and or resolving the fine structure of the flow. The impor tance of the 



meridional flow was demonstrated in protostars and FU Ori stars already by iKley &: Lin 
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( 119961 . Il999l ). who show how matter spread to the poles to completely engulf the accreting 
star. However, in accretion onto a compact star the problem has not been solved because of 
the much smaller scales. 

As the numerical two-dimensional simulations failed to follow the evolution of the ac- 
cretion onto the surface of the compact star, an analytical approach was developed by 
Inogamov fc Sunyeavl (119991 ). who implemented the one-dimensional BL equations with an 
analytic treatment of the meridional direction (assuming a shallow- water equation), where 
the vertical velocity in the boundary layer is permitted to be non-zero. This analytical treat- 
ment was carried out for an accreting neutron-star (NS), an d it was shown that the B L could 
"spread" and cover a significant part of the NS. Later on, (jPiro fc BildstenI l2004al ) applied 
Inogamov & Sunyaev's treatment of the "spread layer" to accreting white dwarf in outburst 
(the optically thick case). The analytical work of iPiro fc BildstenI ( l2004al ) indicates that the 
numerical setup must include sufficient resolution in the radial and meridional planes to cap- 
ture the dynamics and fine structure of the BL. It is the purpose of the present simulations 
to actually provide such a numerical study, by following the evolution of the accretion onto 
(and into!) the surface of the WD in the inner region of the disk and close to the equatorial 
plane (i < 30deg). 

In the one dimensional picture the energy kept by the corotating accreted material per 
unit time is given by eq.4. In two dimension, however, as the matter possibly spreads evenly 
onto the WD surface, its moment of inertia is rather similar to that of a spherical shell and 
the energy kept by the accreted material per unit of time becomes 



1 9 
-\l-mr^ 

2 *3 * 



(6) 



and the remaining energy (eq.4)-(eq.6) 



disk 



(7) 



is available to further spin up the WD. As the matter moves toward the poles, to spread 
evenly on the WD surface, it has excess of angular velocity/momentum (due to the differential 
velocity) which is added to the accreting WD. We therefore see that even for the most 
simplistic two-dimensional model of the boundary layer, the spreading of matter on the WD 
surface involves a significant fraction of the accretion energy (e.g. up to 10% of the disk 
luminosity for a star rotating at l,000km/sec — eq.7). In the case of a star rotating near 
break-up velocity, accretion at the equator adds angular momentum to the star as the matter 
spreads toward higher latitudes. This contributes to a balan cing effect to the spin-dow n of 
the an accreting star rotating near break-up as described by lPopham fc Narayanl (Il99ll ). 
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1.4. The Importance of the Boundary Layer 



The structure of the BL is important to help disentangle and physically characterize 
the different emitting components in an accreting WD system. The details of the emitted 
spectrum depends sensitively on the detailed struc ture of the boundary l ayer and how it 
changes in response to enhanced rates of accretion ( Fisker fc Balsaralboosl ). The boundary 
layer is also important because it may play a role in the still uncertain mechanism which 
drives the outflowing bi-polar winds seen in dwarf novae during outburst and in nova-like 
variables during their high optical brightness states. Our boundary layer simulations will 
eventually provide the theoretical framework required to physically interpet the X-ray EUV 
and FUV emission observed in compact binaries containing accreting degenerate stars. 

However, the structure of the BL is also important because it determines how the ac- 
creted matter ultimately distributes itself in the envelope of the WD, which is important for 
classical novae. In classical novae, the thermonuclear runaway (TNR) ejects part or all of 
the envelope. The initial CNO composition of the burning material should be strongly en- 
hanced compared to t he accreted material to account for composition of the observed ejecta 



(IStarrfie 



see 



deta 



I972I ). Several mechanisms have been suggested for this CNO enhancement 



Josg (120051 ) and references therein). They can be roughly divided into pre-burst mixing 



between accreted material and the underlying CO rich W D (IKippenhahn fc Thomas! Il978 



MacDonaldlll983l : iRosner et al.ll200ll : lAlexakis et al.ll2004l ) by accretion driving instabilities 
or mixing with the underlying material when the conve ctive zone of the thermonuclear run- 



away extends deep enough to dredge up CO material (IStarrfleld et al.l Il972l : iGlasner et al. 
1997h . 



However, before an exploration of the radiation emission characteristics of the bound- 
ary layer can be carried out with full radiation hydrodynamics, we must understand the 
dynamical processes that lead to the formation of the boundary layer. In this flrst stage our 
goal is to calculate the dynamical structure of the boundary layer with suflicient numerical 
resolution to capture the dynamical evo l ution of the accretion flow and its interaction with 
the stellar surface (see lFisker fc Balsaral (120051 )). and we leave the radiation- hydrodynamical 
study for later. 

Our model is presented in section [21 the results are given and discussed in section |3l 
Section m presents the conclusions. The equations we are solving are written down explicitly 
in Appendix |X] while the initial and boundary conditions are described in Appendix [Bl 
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2. Computational model 

The source of the angular momentum transport in the BL probably involves a com- 
bination of magnetic fields and turbulence. Since the angular velocity in the BL is not a 
decreasing functio n of radial distance from the star, it is not clear that the magnetorotational 
instab ility (MRU (IBalbus et al.lll994r) provides angular momentum transport in the BL. Ac- 
tually, Brandenburg et al. ( 19961 ) found that the effective alpha viscosity parameter a would 
go to zero for a rotation law Q ~ r^'^ when g < 0. Nor is the MRI essential at the BL because 
the accretion disk can directly exchange angular momentum and mass with the outer layers 
of the star if an efficient coupling mechanism is found between the disk and the star. Th e 
source of viscosity at th e BL is still debated in th e literature fseelPopharn fc NaravanI (119951 ) ; 
Piro &: BildstenI ( 2004al ): iGodon &: SionI ( 20051 )). Ilnogamov &: Sunyeav Jl999h assumed that 
the viscosity is due to purely hydrodynamical turbulence as in the deceleration of subsonic 
or supersonic flow above a solid surface. The case f or a purely hydrodynamical t urbulence 

however 
e.g. 



DubruUe 1993) 



Orszag fc Kell 



19801 ) 



in the boundary layer was already debated earlier (iZahn 1 119901 : 
the physical process is likely to take place in three dimensions 
by means of streamwise vortices (e.g. [Hamilton fc Habernathy I ( 
instability in a flow over a curved surface, is that of the boundary layer flow over a concave 
surface unstable to the centrifugal instability (as it violates Rayleigh's criterion for stabil- 



19941 )). The most common 



ity — (IRayleigh 



19161)). Th is instability leads to turbulence through the formation of the 



Gortler vortices (jSaric Ill994l ). which tap energy from the laminar flow and poor it into the 
turbulence (this is a non-linear instability leading to a subcritical transition to turbulence). 
However, the boundary layer flow on a convex surface is not subject to this instability, rather 
the opposite, the centrifugal force in such a flow is "restoring". Therefore, the instability 
that is the most likely to take place in the star-disk boundary layer is the Kelvin-Helmhotz 
shearing instabi l ity w hich will occur for a sufficiently large shear (Richardson criteria, e.g. 
Drazin fc Re"id1 Jl98lh l 



Following (jShakura &: Sunyaevlll973l ) , we parametrize the efficiency of the angular mo- 



mentum transport with an a viscosity prescription. In a multidimensional calculation a has 



(Papaloizou &; Stanlev 


1986 


) and ( 


Kley 


1991) 



used z/ = aCsmin{H, hp) where H is the scale height of the disk and hp is the local pres- 
sure scale height in the boundary layer. For these calculations, the pressure scale height 
in the boundary layer that develops on the surface of the star is always smaller than the 
scale height of the disk. The compressible Navier-Stokes equations themselves are written 
explicitly in the Appendix |Al Describing the angular momentum transport with a simple 
shear coefficient means that the dy namics follows the Navier- Stokes equations. Here the 



Navier-Stokes equations as given by iMihalas &: MihalasI (119841 ) are solved in spherical co 
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ordinates {r,6,(j)) on an axisymmetric mesh with 384 ratioed zones in the radial direction 
and 128 ra tioed zones in the merid ional range spanning to 30 degrees from the disk plane 
- same as iFisker fc Balsaral (120051 ) . The star was taken to be a non-rotating O.GMq WD 
with a radius of 9 x lO^cm . The radial extent of our computational domain extended from 
8.9 X lO^cm to 1.1 X lO^cm. The r-directional zones were concentrated towards the inner 
radial boundary with each zone being 1% larger than the previous one. The ^-directional 
zones were also ratioed with the smallest zones being closest to the equator and each zone 
being 1.9% larger than the previous one. Such a ratioed zoning makes it possible for us to 
capture five scale heights of the star's atmosphere as well as the vertical structure of the disk . 



The zone ratioing also concentrates zones around the disk-star interface. Unlike (lKleylll99ll ) 
, who used a mesh with 85 zones in each of the radial and meridional directions, our mesh 
has substantially better resolution. The scal e heights o f WD atmospheres ar e now known 



to be substantially smaller than estimated in (jKleylll99ll : iKley &: Henslerlll987l ). making the 



larger meshes used in this work more essential. The higher resolutions, as well as the use 
of higher order Godunov methods, enable us to substantially reduce the role of numerical 
diffusion. 

In this paper, we wish to study not just the spin-down of the accretion disk but also 
the spin- up of the star. For that reason, we retained five scale heights of the WD's outer 
atmosphere. Care was taken to resolve each stellar scale height in the radial direction with 
at least twenty zones, which ensures a numerically accurate, well-resolved and stable stellar 
atmosphere. To ensure good resolution of the disk's structure, we retained at least fifty zones 
across a disk scale height in the meridional direction. The physical conditions describing the 
structure of the disk are given in Appendix [Bl Typical surface temperatures for hydrody- 
namically accreting WDs are ~ 30,000K and typical disk temperatures are usually taken to 
be ~ 100,000K. This choice of temperatures would make the scale heights of the disk and WD 
atmosphere too small to be resolved with the above-mentioned resolutions. For this reason, 
we systematically allow the temperatures of both the disk and WD atmosphere to be one 
order of magnitude larger than the physical values. As a result, the simulations were carried 
out with a stellar temperature of T=300,000K and a disk temperature of T=1,000,000K. 
While this might seem be a hot temperature for a disk in a CV, it is still much less than the 
virial temperature and the corresponding disk thickness is H/r = 0.03 instead of H/r = 0.01 
for a lO^K disk. The sound speed used in the formulae for the viscosity u was derived for 
the temperatures used in the simulations. We assume a very tenuous and very hot halo 
which consists of the same material as the disk (and the star). The mass in this halo is 
extremely small, making the halo dynamically unimportant and its sole use is to provide 
pressure balance at the upper boundaries of the accretion disk and the star. The fluids that 
make up the disk, halo and stellar atmosphere were tagged with passive scalars and we are. 
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therefore, in a position to assert that the halo gas simply provides pressure support with 
minimal mixing into the gas that makes up the disk or the stellar atmosphere. Such a model 
wi th a disk and h alo that are in dynamical equilibrium with each other was briefly described 
in iBalsaral (120041 ) . Since that previous work did not include the stellar atmosphere, appendix 



IB] describes the physical conditions that were assumed to set up the WD atmosphere, the 
accretion disk and the halo. Symmetry across the disk plane is assumed. 



Observations, and their interpretation in the context of the iPopham fc NarayanI (119951 ) 
model have suggested th at the interesting valu es of a range from 0.1 during outburst to 
0.001 during quiescence. iGodon fc SionI (120051 ) derived a = 0. 004 in the boundary layer 



2003f ). For 
0.01, a = 0.005, 



from the XMM-Newton observations of VW Hyi in quiescence (iPandel et al. 
that reason, we performed a parameter study with a = 0.1, a = 0.03, a 
and a = 0.001 . The same value of a was used in the boundary layer and in the disk. This 
would seem to be a very reasonable assumption if the sam e physical mechanism ( such as a 
thermal instability, ICannizzd (119981 ) or magnetic instability iLivio fc Pringld (119981 )) were to 
make the disk fluid and the boundary layer fluid turbulent. Our simulations do not include 
the effect of radiative transfer in the boundary layer. It is worth noting that in the optically 
thick limit, the BL retains most of the energy that is generated by the viscous stresses. As 
a result, by retaining the viscous energy generation terms in eq. [A5] and assuming that the 
alpha viscosity parameter is large [a ~ 0.1) we mimic the situation where the boundary layer 
is optically thick during outburst. Simulations were also carried out by excluding the viscous 
energy generation terms in eq. |A5j and assuming a ~ 0.001; in that situation we mimic the 
limit where the boundary layer is optically thin during quiescence, i.e. the BL radiates away 
all the heat that is generated by the viscous stresses. Our results, therefore, bracket the two 
extreme cases. Should the high a simulations produce BLs with optical depths in excess 
of 50 or 100 , the optically thick runs would flnd direct applicability to the astrophysical 
problem. Likewise, should the low a simulations produce optically thin boundary layers, the 
optically thin runs would flnd direct applicability to the astrophysical problem. We show in 
the course of this work that such trends are indeed observed in the simulations. In future 
work, we will include a treatment of the radiative transfer terms in the BL thus obtaining 
results that are free of current limitations. The full range of simulations that we report on 
here with the various values of a and the inclusion or exclusion of viscous energy generation 
terms are listed in Table 1. 

Since the inner radial boundary extends into the star, we enforced reflective boundary 
conditions at that boundary. While such a boundary condition would reflect any waves that 
reached the star's surface, we flnd that in practice the waves do not propagate to this depth 
in the present simulations. The use of such a reflective boundary condition represents a 
compromise. While it might reflect back waves that propagate more than flve scale heights 
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into the star, our experience has shown that surface waves almost never propagate to that 
depth. The present boundary conditions do have the positive consequence that they prevent 
any unexpected mass or momentum inflow into the computational domain from the rest 
of the star. The outer radial (computational) boundary was designed to respond to inflow 
or outflow of fluid at the outer open boundary, and it was therefore trea ted by imposin g 
the boundary conditions on the characteristics of the flow as described in (iGodon Ill996bl ). 
Thus outer boundary conditions were imposed directly on the incoming characteristics, and 
computed values from the variables inside the domain were imposed (propagated) on the 
outgoing characteristics. The strategy works well, especially when combined with Riemann 
solvers which also work on the same principle of resolving the inflowing and outflowing 
waves. The equatorial boundary condition in the 6'-direction was reflective (symmetric) and 
the boundary condition at the other 6'-directional boundary was specified as continue. 
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model 


a 


$ 


vl 


0.1 


yes 


v2 


0.03 


yes 


v3 


0.01 


yes 


v4 


0.005 


yes 


v5 


0.001 


yes 


nvl 


0.1 


no 


nv2 


0.03 


no 


nv3 


0.01 


no 


nv4 


0.005 


no 


nv5 


0.001 


no 



Table 1: This table shows our ten models. The first column indicates the model name. The 
second column indicates the a-viscosity employed in the model and the last column shows 
whether dissipative heating was included in the model. 
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m 



The spatially and t e mporally second or der algorithms in RIEMAN N have been descr ibed 
Roe fc Balsaral Jl996h: ISalsaral Jl998aljbl): ISalsara fc Spiceil JlQQQ^ j bll: ISalsaral tooi) and 



Jiang fc Shul Jl996h : ISalsara fc Shu 



usemany ideas from higher order WENO schemes (see 

(I2OOOI )) to reduce dissipation. The matter in the model is subject to the central gravitational 
field of the underlying WD. The model uses an ideal gas (7 = 5/3) of a fully ionized solar 
composition (/i = 0.62g/mole) and assumes no radiative transport. 



3. Results and Discussions 

In the next subsections, we focus on several important aspects of the boundary layer 
dynamics as follows. In section 13.11 we focus on the density profile of the boundary layer. In 
section 13.21 we consider the pressure profiles in the boundary layer. In section 13. 3[ we focus 
mainly on the the evolution of angular momentum in the boundary layer . In section 13.41 
we check the stability of the flow in the boundary layer. In section 13. 5[ we discuss accretion 
based instabilities and in section 13.61 we relate the computations to observations. 



3.1. Density Structure of the Boundary Layer 

Figs la and lb show the logarithm of the density in the boundary layer with a = 0.1 
and a = 0.001 respectively. Figs Ic and Id and do the same for the runs with the same 
values of a but with the viscous heating switched off in eq. |A5j . The full extent of the 
computational domain is shown in the 6'-direction. Only a small portion of the inner radial 
direction is shown and the values on the x-axis of the plot are in units of lO^cm, making it 
possible to measure the radial coordinate . The same convention for labeling figures applies 
to all other figures in this paper where fiow variables are imaged. 

In all cases, the poloidal velocity is overlaid as vectors, enabling us to trace the fiow 
of matter on the surface of the star. Thus Fig. la pertains to the optically thick limit 



Fig. 1. — (a) left panel. The logarithmic density for model vl is shown after one Keplerian 
rotation, (b) right panel. The logarithmic density for model v5 is shown after one Keplerian 
rotation, (c) next page - left panel. The logarithmic density for model nvl is shown after 
one Keplerian rotation, (d) next page - right panel. The logarithmic density for model nv5 
is shown after one Keplerian rotation. The full extent of the computational domain is shown 
in the ^-direction. Only a portion of the radial direction is shown and the values on the 
X-axis of the plot are in units of lO^cm , making it possible to measure the radial coordinate. 
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(outburst state) while Fig. Id corresponds to the optically thin limit (quiescent state). We 
see that a thick, dense boundary layer forms in both Figs, la and Ic while that is not the 
case in Figs, lb and Id. Thus, physically thick boundary layers form in outburst and the 
result is independent of whether the viscous energy generation is included in the energy 
equation. The boundary layers that form in quiescence tend to be physically thin. This 
shows that the viscosity is the primary discriminant in determining the structure of the BL. 
The poloidal velocity vectors in Figs la and Ic also show us that the velocity increases with 
increasing distance from the star's surface. The velocity vectors in Figs, lb and Id show 
a similar trend though it is harder to see because of the smaller physical extent of the BL. 
Such a velocity structure is also known to occur in terrestrial fluid dynamics when a viscous 
fluid flows over a stationary solid surface, forming a boundary layer at the solid's surface. 
This shows us that our simulations are producing a valid result which is consistent with our 
intuition. It also shows us that our decision to refer to these star-disk layers as boundary 
layers is well-motivated. 

The accretion rate is high enough in the a = 0.1 cases that the infalling matter depresses 
the stellar atmosphere close to the equator. While this is not so evident in Fig. la-Id, we 
will show in section 13.51 that the infalling matter can excite gravity waves on the surface of 
the star. Likewise, in section 13.31 we will show that the high a case spins up a significant 
portion of the star's atmosphere while itself being spun down. The net result of this process 
is that the toroidal velocity of the disk becomes sub-Keplerian at larger radii from the star 
as a increases, resulting in broader boundary layers. Fig. 2 shows the mass accreted as 
a function of time. This figure was generated by integrating the mass from the disk-star 
interface to the top of the BL where dVfp/dr = 0. Fig. 2 shows us that the accretion rate is 
higher in the high-a cases, as expected. Fig. 2 also shows us that for the high a cases the 
accretion rate seems to drop off with time. This is entirely a result of the fact that this round 
of simulations does not include radiative transfer. As a result, the base of the boundary layer 
heats up, with a corresponding increase in pressure. The radiative cooling times in a real 
accretion disk are short enough (i.e. even smaller than an orbital time) that the boundary 
layer would cool down quite rapidly. The consequent pressure reduction would then permit 
accretion to proceed unimpeded. 

Because of resolution constraints, we had to use temperatures that are somewhat larger 
than those that prevail on accreting white dwarfs. In subsequent work, we intend to overcome 
this constraint. It is, nevertheless, interesting to relate physical depth to optical depth of 
the boundary layer. In doing that, it is important to exclude the disk fluid that is spinning 

Fig. 2. — The mass accretion is shown as a function of time for models vl-v5. 
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close to the Keplerian velocity. We, therefore, define the disk-boundary layer interface as 
the region where the gradient of the angular velocity vanishes: dVt/dr = at r = r^. In the 
boundary layer (r < r^) the matter is spinning with a toroidal velocity that is smaller than 
the Keplerian velocity. For that sub-Keplerian accreted material, we plot the optical depth 
due to Thompson scattering as a function of meridional angle in Fig. 3. Here we define the 
optical depth as — updiskdr, where k = 0.34 gcm~^ is the Thomson scattering opacity of 
a fully ionized solar composition. This constitutes the minimum amount of scattering and 
thus provides a lower bound of the opaqueness of the matter. 

The boundary layer is technically defined by the radius where the radial gradient of the 
accretion fiow's angular velocity disappears. We use that definition for the rest of this paper. 
Because we track the fiuids that make the disk, halo and star, we are in a position to isolate 
just the disk material that is within the boundary layer. The optical depth of this material 
(in the radial direction) is very important because it sets the emission characteristics of 
the accretion belts that have been identified in the observational literature. Fig. 3 shows 
us that physically thick boundary layers in the high-a limit also tend to be optically thick 
while physically thin boundary layers in the low-a limit tend to be optically thin. The 
physical implication of that is that during outburst the BL could suppress the emission from 
a fraction of the star's surface. If the BL is optically thick with an optical depth ^ 1 then 
the hottest part of the boundary layer, which prevails at the disk-star interface would not be 
visible. However, during both the transitions, from outburst to quiescence and vice-versa, 
it is possible that this hot layer might become visible, giving one a direct view of boundary 
layer heating. Th is might show up as an additional UV component such as the one seen by 



Sion et al.l (120051 ) in VW Hydri during its transiti on from qui e scence to outburst. A simila r 



UV component has been detected in U Gem by iLong et al.l (119931 ) ; iFroning et al.l (120011 ) ; 



Szkody et al.l (120021 ) during its transition from outburst to quiescence. Obtaining matched 
measurements of the velocity and UV excesses would allow one to further corroborate the 
scenario presented here. We see that the model in Fig. la produces an optically thick 
boundary layer during outburst. It is, therefore, consistent with our claim in Section 2 that 
inclusion of the viscous dissipation term in eq. [A5] provides a rather realistic representation 
of the structure of the boundary layer in outburst. Likewise, the model in Fig. Id, results 
in an optically thin boundary layer during quiescence. It is thus consistent with our claim 
that excluding the viscous dissipation term in eq. [A5] is similarly more representative of the 
structure of the boundary layer in quiescence. 



Fig. 3. — The optical depth along the radial direction is shown for models vl-v5 as a function 
of the angle from the equatorial plane. 
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3.2. Pressure in the Boundary Layer 



Figs 4a and 4b sliow tlie logaritlim of tlie pressure in the boundary layer with a = 0.1 
and a = 0.001 respectively. Figs 4c and 4d do the same for the runs with the same values 
of a but with the viscous heating switched off in eq. |A5] . 

From Figs. 4a and 4b we see that the pressure of the accreted fluid is highest in the 
equatorial plane. This high pressure can be attributed to a combination of viscous dissipation 
as well as the ram pressure due to infall of accreting material. The pressure gradient in the 
meridional direction on the surface of the star, therefore, drives the meridional flow that 
develops on the WD's surface. The models that were used in Figs. 4c and 4d did not include 
viscous heating in eq. |A5j . They, nevertheless, show that the pressure of the accreted fluid 
is highest in the equatorial plane and in this instance, the increased pressure is entirely 
attributable to the ram pressure due to infall of accreting material. We also see that the 
equatorial pressure increase in Fig. 4b is less than that in Fig. 4a and, similarly, for Figs. 4b 
and 4d. The smaller pressures in Figs. 4c and 4d relative to Figs. 4a and 4b can be attributed 
to the additional contribution from viscous heating. Fig. 4a-4d, therefore, serves to show 
us that the formation of the pressure gradient in the meridional direction, which then drives 
the meridional flow in the boundary layer, is a very commonplace phenomenon. In other 
words, any flow that is put in a similar geometry and is made to experience similar viscous 
stresses would naturally form a similar boundary layer, showing the ubiquity of boundary 
layer formation. (The only other requirement that such a BL flow has to satisfy as the disk 
material migrates polewards on the star's surface is the ability to lose angular momentum. 
The next sub-section will show that it can accomplish this very effic iently by spinning up 
the un d erlying layers of t he stel lar atmosphere.) Thus the decision by llnogamov fc Sunyeav 
( 119991 ): iPiro fc Bildsteru (j2004al ) to include a non-zero meridional velocity in their models 
was of central importance in forming the inherently multi-dimensional boundary layers that 
we see in our simulations. Such boundary layers also form in ac creting neutron stars, an d 
TTauri stars that are go ing through the FU Orionis phenomenon iKley fc Lin I (119961 . Il999l ) ; 



Balsara fc Fisker I (j2005|). 



Fig. 4. — (a) left panel. This figure shows the logarithmic pressure for model vl after one 
Keplerian rotation, (b) right panel. This figure shows the logarithmic pressure for model v5 
after one Keplerian rotation, (c) next page - left panel. This figure shows the logarithmic 
pressure for model nvl after one Keplerian rotation, (d) next page - right panel. This figure 
shows the logarithmic pressure for model nv5 after one Keplerian rotation. 
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3.3. Toroidal and Poloidal Velocities 

As the simulations start, the shear force transfers angular momentum between differen- 
tially rotating disk annuli so that matter can move inwards and accrete on the star at the 
footpoint of the disk. Angular momentum transfer between the disk material and the stellar 
surface is also necessary so matter can move to higher latitudes. Otherwise, the centrifugal 
barrier prevents matter from leaving the footpoint of the disk. It is this transfer of angular 
momentum which spins up the existing surface layers of the star. Angular momentum is also 
directly advected onto the star, because the orbiting disk material eventually accretes and 
forms the new surface. 

Fig.Sa shows the resulting specific angular momentum after 3/4 of a Keplerian evolution 
for the a = 0.1 case, whereas fig.Sb shows it for the a = 0.001 case. Fig.Sa shows that the 
halo is spun up as is the underlying star at the footpoint, where the disk connects with the 
star. Moreover, we observe that by this time in the simulation the disk material has spun 
up all five atmospheric scale heights of the star at equatorial latitudes. At higher latitudes 
the spreading disk material is less dense and thus takes longer amounts of time to drag the 
star's surface around with it. In fig. 5b, which is at the same time as fig.5a and uses the same 
color scale, we see that the disk has not spun up the star's atmosphere even at the stellar 
equator. 

Fig.6 shows the toroidal velocity and sound speed in the disk's midplane as a function of 
radius for the simulated alpha-viscosities. The cases of a = 0.01, a = 0.005, and a = 0.001 
do not show any significant spin-up of the star after 3/4 of a Keplerian rotation, whence there 
is not significant shear connection between the disk and the stellar surface. Therefore disk 
matter retains its Keplerian motion much closer to the star which means that even though 
the accretion rates are smaller, significant amounts of the dissipation can still be generated 
in the BL close to the star. 

The matter in the layers of the star that we simulate is indeed non-degenerate. Even so, 
the high temperatures cause it to have a rather low specific molecular viscosity. As a result, 
turbulence and/or threaded magnetic fields might nevertheless be the two most prominent 
ways to transfe r a significant amount of angular momentum to the deeper layers of the WD 
JPurisenlEgTsh . where the relative motion drives the mixing between the accreted surface 



material and the deeper layers. A supersonic component in the toroidal direction over most of 



Fig. 5. — (a) left panel. A color plot of the specific angular momentum, Q = v^/r, for the 
a = 0.1 case after 3/4 of a Keplerian rotation period, (b) right panel. A color plot of the 
specific angular momentum for the a = 0.001 case after 3/4 of a Keplerian rotation period. 



- 20 - 



the BL is obtained for all values of a (see fig. 6). The supersonic toroidal velocities mean that 
the flows are susceptible to the rapid development of gravity wave and/or Kelvin-Helmholtz 
instabilities in three dimensions. Such instabilities can be studied by taking a small slice 
of the interface from our simulations and extending it in the toroidal direction in a three 
dimensional simulation. Our present simulations could then provide the velocities which 
drive such instabilities. The calculation of the turbulent mi xing which obtains from these 



instabilities must therefore be calculated using other mod els ( iKippenhahn fc Thomas! 1 1978 
MacDonaldlll983l : iRosner et allboOll : klexakis et aDl2004l ). 



Fig. 7a shows the poloidal Mach number of the accreted material on the surface of the 
star for a = 0.1. Fig. 7b shows the same for a = 0.001. The same scale is used for both 
figures. We see that the Mach number is mildly transonic in Fig. 7a while it is subsonic 
in Fig. 7b. Only the a = 0.1 case has a transonic component in the poloidal directions, 
whereas the poloidal f lows for a = 0.03, a = 0.01, a = 0.005, and a = 0.001 remain subsonic 



see figs. 5a and 5b) (IFisker fc Balsarall2005l ). 



3.4. Hydro dynamic Instability in the Boundary Layer 

The present simulations are based on an alpha-viscosity formulation, which implicitly 
assumes an underlying model for the sub-scale turbulence. However, the material that 
accretes on to a white dwarf has very low viscosity. Observations have not revealed the 
existence of magnetic fields at the WD surface. Besides, the radial gradient of the angular 
velocity has the wrong sign for the MRI to act. It is, therefore, interesting to explore the 
role of other hydrodynamical instabilities at the surface of the star. It is important to 
remember that the alpha-viscosity formulation does smear the velocity gradients. Yet, if 
these gradients are amenable to the development of a persistent hydrodynamical instability 
then the simulation should reveal its existence to us. It is in that spirit that we try to identify 
an instability that might give rise to a sub-scale hydrodynamical turbulence on an accreting 
white dwarf's surface. 

There are several instabilities that can possibly lead to turbulence in boundary layers, 
such as the centrifugal instability (Rayleigh's criterion), the shearing Kelvin-Helmholtz in- 
stability, or even the Rayleigh- Taylor instability (buoyancy forces). Depending on the exact 



Fig. 6. — Toroidal (0) velocity in the disk's midplane shows the extent of the boundary layer 
as well as the sound speed after 3/4 of a Keplerian rotation period. Note that asymptotes 
to a higher value than Cg^a- 
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angular velocity profile, pressure and density gradients each of the forces involved compete, 
some are stabilizing while other are destabilizing. The sufficient condition for linear stability 
of a rotating , radially stratified fiuid under the infiuence of (a radial) gravity was given by 



Sung I fll974[ ) 



1 ^ kf lm2 fdnY 

where, m and kz are the angular and vertical mode numbers (in cylindrical coordinates 
{zj, 0, 2;)), M is defined as M = /c^ + rr? jw'^ ^ geff is the effective gravity defined as 



9eff 



GM 



S'ro is the Schwarzschild discriminant (ISchwarzschild Ill906l ) 
in which 



dp 
dw 



ad 



dp 
dw 



dp 
dw 



ad 



p dp 
Vip dw 



is the adiabatic density gradient; and $ is the Rayleigh discriminant (IRayleigh III88OI . Il916l ) 



1 d /^2q\2 

w'^ dw 



In the present case we do not consider vertical modes kz and only consider the equatorial 
fiow for which r = w. There are two distinct cases as follows. 



The axi-symmetric case (mode m = 0) gives the Solberg-H0iland criterion (ISung Ill974 



19751) 



1 



geffSr + $ > 0. 



The perturbatio ns of an accretion disk to a convective instability has already been stu died 
in this context (ILivio &: Shaviv 1119771 : iRiidiger et al. 1120021 : 1 Johnson &: Gammie II2006I ). In 
the limit of vanishing rotational velocit y, the Solberg-H0iland criterion simply leads to the 
classical Schwarzschild condition Sr > (ILebovitz Ill965l . 119661 ) , which is also the condition in 



Fig. 7. — (a) Left panel. This figure shows the poloidal Mach number with poloidal velocity 
vectors overlaid for model vl after one Keplerian rotation, (b) Right panel. This figure 
shows the poloidal Mach number with poloidal velocity vectors overlaid for model v5 after 
one Keplerian rotation. 
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the vertical direction in the disk (iLivio fc Shaviv 1119771 ). In the hmit of a vanishing gradient 
of the density, the Solberg-H0iland criterion simphfies to Rayleigh's criterion for rotating 
fluids. The axi-symmetric instability in the boundary layer could possibly lead to a ring-like 
structure oscillating and in many ways similar to the "breathing" mode of an unstable star. 
This mode will not lead to turbulence, but could it explain the short period oscillations 
observed in CVs? The mode would depend on the density, temperature and rotation rate 
(Solberg-H0iland criterion), and its period and coherence would vary with these parameters. 
However, dwarf nova oscillations (DNOs) in CVs exhibit a ISOdeg pha se jump at eclips e 
( iPatterson Ill979l ). and frequency doubling (presence of the first harmonic [Patterson I (jl98ll ). 
and therefore they cannot be explained by an m = mode. On the other hand. Quasi- 
periodic oscillations (QPOs) in CVs do not all exhibit the same characteristics (some are 
believed to form at larger radii in the disk), and, therefore, we cannot completely rule out 
the m = mode in the boundary layer as an additional mechanism to produce quasi-periodic 
oscillations (QPOs). 

We now turn our attention to the non-axisym metric case. In the case m 7^ 0, a sufficient 
condition for stability that does not involve kz is (ISung Ill974l ) 



1 

-■i 

4 



p " 4 \dr J 
This condition is usually written in the form of a modified Richardson number 

I L ap L ap \ I a\L \ ^ 

9eff 



(9) 



Ri 



1 dp 1 dp 
Tip dr p dr 



dr 



> 



1 



(10) 



Thi s is a gener a lizati on of the Miles-Howard theorem (IMiles I (Il96ll ): iHoward I (Il96ll ): see 
also Chinionas (Il970 ) ). Th e Miles-Howard theorem itself reduces to the original Richardson 
criterion (IRichardson Ill920l ) when compressibility is omitted and the Schwarzschild discrim- 
inant is replaced by the buoyancy term alone. The modified Richardson number was consid- 
ered in a few ana lytical studies to analyze the stability of material accreting on the surface 
of a white dwarf (jDurisenlll977l : iKippenhahn fc Thomaslll978l : iLivio fc Truran 1119871 ). Here, 
we have the opportunity to use results from numerical simulations to evaluate the modified 
Richardson number in the b oundary layer. Using the definition of the Ledoux discriminant 
( ILedoux fc Walraven Ill958l ). 

^ 1 d In P d In p 
Fi dr dr ' 

the (modified) Richardson condition for stability is then written 



Nlv 



rfi^y^ 1 



(11) 
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wher e A is related to the buoyancy (or Brunt-Vaisala) frequency Nbv f Pesnell Ill986 : Livio fc Truran 
H, e.g.) 

9effA. 



Ar2 



The the flow is unstable for i?j < 1/4. In Fig. 8 we show the Richardson number, Ri in the 
0-direction (eqs. 10, 11), for run 91 and 95. The Richardson number is smaller than 1/4 in 
most of the domain and therefore the flow is unstable. It is i mportant to remark that while 
the shear could stabilize the f low (Johnson fc Gammie internal gravity waves can be 

reflected from a shear laver (IVan Duin fc Kalder Ill982h . and can be over-reflected from a 
rigid boundary ( Shachdev fc Satya Narayanan I 19821 ). It has also been shown that modes 
can be unstable at the star-disk interface due to the propagation through the corotation 
(ITsang fc Lai 1120091 ). Unstable modes could be trapped and grow between (i) the surface of 
the WD and the strong shear region; and/or (ii) between the strong shear region and the 
corotation radius at larger radii in the inner disk, in a ma nner similar to the Papaloizou- 
Pringle instability (jPapaloizou fc Pringle Ill984l . Il985l . 119871 ) observed in simulations of disks 
with a rigid inner boundary (iGodon Ill997al ). The question of how the instability will develop 
cannot be addressed without full 3D simulations, which are beyond the scope of this work. 
However, a Richardson number < 1/4 in the boundary layer region and at the stellar surface 
raises the possibility of an instabili ty. This instabil i ty coul d develop in the form of waves 
in the "spread layer" as studied by iPiro fc BildstenI (l2004bl ) to explain DNOs in CVs, and 
raises the possibility of strong mixing between the hydrogen rich freshly accreted material 
and the WD material. 



3.5. Formation of surface gravity vi^aves 

For the a = 0.1 case, the outer edge of the boundary layer (where dQ/dr ~ 0) is 
located at l.OGr^,. At r = 1.03r=K, is about 92% of vk- This means that only 84% of 
the gravitational pull is supported by the centrifugal force while the rest is supported by 
pressure. The pressure for a = 0.1 is illustrated in fig.4a and comes from the build up 
of density due to the high accretion rate facilitated by the high viscosity which drags down 
material from the disk and also keeps it from moving rapidly towards the poles once it makes 
contact with the WD surface. This results in a dense band at the foot point of the disk which 
causes a surface gravity wave of surface matter to spread towards the poles. Such surface 



Fig. 8. — The toroidal Richardson Number as a function of the radius r in the equatorial 
plane {6 = 90°), for run 91 and 95. The stellar radius is located at r = 9 x 10*^ cm. In most 
of the domain the Richardson number is smaller than 1/4. 



-24- 



gravity waves have been studied in the context of terrestrial phvsi cs bv iMiles I (119571) and 
the role of turbulence in exciting these waves has also been studied IPhilhps I JiOStT and the 



simulations show a similar phenomenon on the surface of strongly accreting white dwarfs. 
The matter inflow is, however, sufficiently large to cause the disk material to overflow the 
gravity wave supersonically as described above and shown in fig. 5a. Our present simulations 
do not show any evidence of wave breaking. However, should wave breaking reveal itself in 
simulations where the accretion is more vigorous, it would provide a further mechanism for 
directly mixing accreted material with the CNO-rich material that lies beneath the surface 
of the white dwarf. 

The propensity for the creation of gravity waves and their magnitude decreases rapidly 
with lower values of a. Even for a = 0.03, the effect is only marginal and for lower values it 
is no longer present. The effect is thus only present for large values of a, and possibly during 
dwarf novae in the outburst stage. Furthermore the gravity wave might be transient as the 
surface adjusts to a fluctuations in the accretion rate. For lower values of a, the matter 
accretes slowly inwards and spreads towards the poles in a uniformly thick layer. 



3.6. Dissipative heating of the BL 



The local rate of ene rgy dissipation depends on the divergence of the angular velocity 
( iMihalas fc Mihalaslll984l ). so flg.6 also indirectly shows where heat is released. Fig. 9a shows 
that the supersonic impact of the accretion flow dissipates the most energy. However, energy 
is also dissipated at the interface between the BL and the surface. This is better seen in 
Fig. 9b which shows that the toroidally rotating matter dissipates rotational kinetic energy 
along its entire interface with the surface as it slides against it. 

In situations where the BL is optically thin, the radiation will be emitted from the very 
same regions where the heat is generated in flg. 6b. This may explain the optically thin 
BL and line emis s ion s een in observations of U Gem during the quiescent phase (see also: 
Fisker &: Balsaral (120051 )). Even in outburst. Fig. 3 shows that the optically thick boundary 
layer has an optical depth of 100. We therefore expect that the radiative diffusion of 
the heat in the radial direction will be faster than its advection in the poloidal direction. 
Consequently, our simulations provide a dynamically consistent reason fo r the forma t ion o f 
a warm belt of accr eted matter, as was s uggested in the observations of ISion et al.l (120051 ) 
and the modeling of iGodon fc SionI (120051 ). 
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4. Summary and Conclusion 

We have numerically simulated the structural dynamics of the BL for an accreting white 
dwarf surrounded by an a-disk for different values of a. The structure and dynamics of the 
BL is important, because it determines the specifics of the radiated spectrum that is emitted 
from this region. The energy radiated from the BL could account for up to half of the total 
energy released in the accretion process. For high values of a, the BL is optically thick and 
extends more than 30 degrees to either side of the disk plane after 3/4 of a Keplerian rotation 
period (tx=19s). 

The simulations also show that high values of a result in a spreading BL which sets 
off gravity waves in the surface matter. The accretion rate can be high enough in the 
high a cases to cause a depression to form on the surface of the star. The accretion flow 
moves supersonically over the depression that is formed, making it susceptible to the rapid 
development of gravity wave and/or Kelvin-Helmholtz instabilities. The low viscosity cases 
also show a spreading BL, but here the accretion flow does not set off gravity waves and it 
is optically thin. 



We argued in the paper by ISion et al.l (120051 ) that an accretion belt might be sustained 



after a long period of perhaps thousands of dwarf nova events such that an equilibrium or 
steady state is established between the Richardson number and the average rate of accretion. 
Our 2D simulations in this paper has not been followed long enough for any steady state 
behavior to be observed. Therefore, it may be premature for us on the basis of the limited 
extent of these first simulations to argue in favor or against the formation of an accretion 
belt. 

If the boundary layer remains optically thick following an episode of high accretion, 
this could explain the "second components" of FUV flux seen in several dwarf novae during 
quiescence. It should also be pointed out that a hot equatorial accretion region could persist 
even after the material has reached co-rotation with the white dwarf. For example, a hot 
(50,000K) region of low rotation is seen in successive HUT spectra of U Gem taken early 
after an outburst and very late after the same outburst. 

The susceptibility of the flow in the boundary layer to undergo a purely hydrodynamic 



Fig. 9. — (a) Left panel. A logarithmic color plot of the dissipation rate for the a = 0.1 
case after 3/4 of a Keplerian rotation period, (b) Right panel. A logarithmic color plot of 
the dissipation rate for the a = 0.001 case after 3/4 of a Keplerian rotation period. Notice 
that the color scale is different for fig. 8a. 
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instability (modified Richardson number < 1/4) was found assuming an a priori alplia viscos- 
ity parametrization consistent witfi tfiat of tlie disk. Tlie unstable flow would be very effective 
in mixing of the accreted material with the outer layer (few scale heights) of the WD enve- 



lope. However, it has been shown (IPapaloizou &: Szuszkiewicz 1 1 1994 iNarayan et al. I Il994 



Kato fc Inagaki 1994 : Godon 19951 ) that the turbulent viscosity law for a non-Keplerian 
disk cannot be represented by the standard alpha viscosity prescription [ a =constant), as 



alpha is proportional to the shear (iGodon Ill995l : lAbramowicz et al. 1119961 ) 



a oc dVl / dr. 

Even more so, in the boundary la yer, the turbulent viscosity is proportion al to the square of 
the turbulent Mach number Jv(l (IShakura fc Sunyaevlll988l : iGodon Ill995l ). which decreases 
quickly for for subsonic turbulence (supersonic turbulence quickly dissipated due to shocks) 



T^BL = asLCsH, 



ubl = — 



Supersonic ra 
ity as shown 


dial infall in the boundary layer also leads to a smaller turbulent viscos- 


bv (] 


^aDaloizou &; Szuszkiewicz 


1994; 


Naravan et al. 


1994; 


Kato & Inaeaki 


1994; 


Godon 


199. 


5), and the formulation of such a viscosity can also be obtained consid- 



ering causality in orie-dirn ensional steady-state models of accretion disk boundary layers 



(jPopham fc NarayanI Il992l ) , namely the radial infall vi scosity cannot be larger than the 



maximum speed of the turbulence. In that case one has (IGodon Ill995l ) 

Sbl 1 



l^BL = asLCsH, 



OCBL 



a- 



H Mf 

where Adr > 1 is the radial infall Mach number, and it is smaller than the turbulent Mach 



number Adt (see also iPopham fc NarayanI (119921 )) 



In these regions close to the stellar surface (and even more so away from the equatorial 
region) , one can use a viscosity of the form 



V oc 



(12) 



where f* is the stellar rotational velocity, v is the velocity of the flow, and z is the distance 
from the surface of the WD. It is clear that such a viscosity would decrease rapidly as z — > 0, 
and would agree with the fact that the size of the turbulent Eddies be limited by the presence 
of t he stellar surface and their velocity proportional to the change of v over a distance z (see 
also iLandau fc Lifshitzl (119871 )). 



Piro and Bildsten have shown that the properties of the boundary ( "spread" ) layer vary 
depending on the value of the viscosity, however, like us, they assume an alpha viscosity in 
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which alpha is not a function of the coordinates, the radial infall velocity or or the shear. It 
is our aim, in future simulations, to carry out simulation of the boundary layer, by assuming 
a viscosity law that is greatly reduced in the boundary layer. For this purpose, one could use 
a combination of viscosity prescriptions suggested in the above mentioned works. Since the 
results presented here as well as Piro & Bildsten's results depend strongly on the value of a 
(which was assumed to be constant), it is clear that a modified alpha viscosity prescription 
can lead to some new and unexpected results. 

The evolutionary simulations presented here, while followed for only a brief interval of 
time, represent the first successful, fully hydrodynamic treatment of the the flow of accreted 
matter into the white dwarf surface layers. All previous attempts to follow accretion hy- 
drodynamically from the boundary layer into the white dwarf surface layers were either not 
computationally viable or the underlying white dwarf was treated as a solid boundary. 

The successful convergence of our dynamical model simulations opens the door to going 
well beyond this first stage and ultimately follow the dynamical evolution over much longer 
time intervals for both the high viscosity and low viscosity cases with the inclusion of radiative 
processes. 

PG wishes to thank Steve Lubow for a discussion on the importance of the adiabatic term 
in the modified Richardson number, and Mario Livio for his kind hospitality at the Space 
Telescope Science Institute. This work is supported by NASA ATFP grant NNX08AG69G to 
Villanova University and the University of Notre Dame. Participation by EMS and PG was 
also supported in part by NSF grant AST08-07892 and NASA ADP grant NNX04GE78G 
to villanova University. 



A. Numerical model 



The source of the angular momentum transport in the BL is a combination of mag- 
netic fields and turbulence. However, an a priori prescr iption, in particular in the BL, is 
a source of disagreement (see iPopham fc NarayanI (119951 ) and references therein) . Instead 
the efficiency of the angular momentum transport can be parametrized with a coefficient, 
a (jShakura fc SunyaevI Il973l ). Describing the angular momentum transport with a sim- 
ple shear coefficient means that the dynamics follows the Navier-Stokes equations. Here 
the Navier-Stokes e quatio ns are solved in spherical coordinates i = {r,6,(f)) as given by 
Mihalas fc Mihalad (Il984j ). The accelerations are defined by the time-derivative of the ve- 
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locities which are given in spherical coordinates. 
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Here are the accelations, p is the pressure, and ( is the coefficient of bulk viscosity, which 
we set to zero to treat the fluid as a Maxwellian fluid. The kin ematic viscosity v = p/ pis 
set by the alpha-disk prescription of IShakura fc SunyaevI (119731 ) . so v = aCsH, where H is 
the disk scale height which should properly correspond to the turbulent turnover scale with 
convective bubbles moving at sound speed c^. 



The viscous dissipation function is given by 
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The dissipation function enters the heat equation as 
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where T is the temperature, = R/{^ — 1) is the specific heat capacity. As mentioned 
in section 2, the heat equations does not contain transport terms, so g = 0. This means 
that dissipated heat can only be lost be pressure work. Conversely, if the dissipation term is 
neglected as is the case in the optically thin limit, the only way to change the temperature 
at any given point is through pressure work. 



B. Computational domain setup 



We set up an accretion disk around the WD. This setup comprises three parts: the 
disk, the stellar atmo sphere and the halo. This improves on the previous simulations of 
Kley fc Henslerl (119871 ) which did not include a stellar atmosphere. For stability reasons it 
is important that the pressure of the three components match at their respective interfaces. 
This is to prevent shocks from dominating the solution. We have therefore carefully designed 
the domain to ensure the partial pressures are identical at the common interface. 

The radius of the star is r = 9000km. We assume that the disk is symmetric in the 
plane and that is it axisymmetric. We use the model of iBalsaral (120041 ). In their notation, 
we have a = 1/2 and 7 = 5/3, so their model effectively reduces to a constant density in the 
disk plane for our case. We assume the the initial halo, atmosphere and disk are isothermal 
and obey the poljd^ropic equation of state 



P = Kp' (Bl) 
^. Setting a temperature for the disk then fully determines K. 



We set po = 1-2 X 10~^gcm 
Above the disk plane the partial pressure of the disk matter follows the prescription of a 
plane parallel atmosphere where the pressure falls off with the e- f olding distance given by th e 



scale height. The density is calculated according to eq. |Blj (see iLandau fc Lifshitzl (119871 ). 



The temperature of the WD is set to = 300000K. The atmospheric scale height is 
given hy H = RT^/gfi, where g = GM/r'^, R is the gas constant and fi is the mean molecular 
weight. The inner boundary of the computational domain is set to r* — 5H^. The disk scale 
height is given by = Cg^/rjg, where Cg = a/7-P/ p is the local sound speed. The outer 
boundary of the computational domain is set to r* + 3i7^. The disk is set up with a Keplerian 
velocity profile. The disk temperature is set to = lOOOOOOK. Pressure matching between 
the disk and the stellar atmosphere is achieved by selecting so that = Pd at the base of 
the disk at r = r*. The pressure pr ofile of the WD atraosphe re follows the standard solution 
of a plane parallel atmosphere (see iLandau fc Lifshitzl (119871 )) 



In these models we do not have a self-consistent disk corona, which would provide 
pressure balance at the upper boundary of the accretion disk. For that reason, we use a 
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very hot, initially static halo that provides pressure balance both to the upper boundary of 
the disk as well as to the WD's atmosphere. The halo's temperature is chosen to be high 
enough that the mass in the halo is less than the mass in the disk or the WD atmosphere 
by a factor of 35, making it dynamically unimportant. The temperature of the halo is set 
so Th = SOTrf. The scale height of the halo is determined by the gravitational component 
in the z-direction. The density of the halo is set to match the pressure of the disk one disk 
scale height above the disk base at r*. 

The careful pressure matching described above ensures that the dynamical instabilities 
at t = arc negligible. Therefore, the computational model will quickly stabilize to a 
dynamical equilibrium. 
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